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Gravitational constant 

p x 1 vector of constant Lagrange multipliers 
Atmospheric reference density 
Maximizing function 
Terminal constraint vector 
Stopping constraint 

Partial derivative evaluation for nominal estimates 


ACKNOWLEDGEMENT 

The author gratefully acknowledges the cooperation and programming 
efforts of Mr. Bill Keenum, Computer Sciences Corporation, in connection 
with the development of the two program decks in the appendix. 


TECHNICAL MEMORANDUM X -646 36 


APPLICATION OF THE STEEPEST ASCENT OPTIMIZATION 
TRAJECTORY PROBLEM METHOD TO A REENTRY 

INTRODUCTION 


Trajectory optimization problems are usually concerned with the task 
of controlling a dynamical system such that a particular mission trajectory is 
accomplished with some measure of performance being extremized (maximum 
or minimum) . Once this measure of performance is selected, it is used with 
the system equations and the initial and terminal boundary conditions to form- 
ulate the optimization problem. The mathematical details of a particular form- 
ulation depends, in general, upon the complexity of the problem. The methods 
resulting from the various formulations can be grouped into two classes: 

( 1) direct and (2) indirect. Most direct methods are based upon the results 
of Kelly [1, 2] and Bryson and Denham [3] known as the method of steepest 
ascent while the indirect methods stem from either ( 1) the calculus of varia- 
tions [4], (2) Pontryagains' maximum principle [5], or (3) dynamic programm- 
ing [6] . To obtain explicit solutions to the optimization problem using the in- 
direct methods, a nonlinear two-point boundary problem must be solved. This 
difficulty is circumvented when the direct method of steepest ascent is used. 

The purpose of this report is to present an application of the above 
mentioned steepest ascent method to an Apollo three-dimensional reentry 
optimization problem. This particular problem has been investigated by 
Colunga [7] using a modified sweep method (MSM) . The MSM is a second- 
order indirect numerical optimization method whereas the steepest ascent 
is of first order. 

THE Dl RECT METHOD OF STEEPEST ASCENT 

Discussion 

Steepest ascent is an iterative procedure in which the nominal or 
beginning values of the control variables must be supplied by the analyst. 
Optimum values for the control variables are determined through a sequence 
of perturbations to the control variables; i.e. , the control variables are 
perturbed by a certain amount and the resulting first-order predictions of 
changes in the payoff and constraint functions are determined. The steepest 
ascent method then seeks the perturbed control variable time history which 



results in maximizing or minimizing the payoff function while simultaneously 
satisfying the constraints. The steepest ascent theory can be summarized 
as follows: If one goes through a sequence of the prescribed cycles, which 
results in improved trajectories, then eventually a trajectory is obtained 
which is reasonably close to the optimum. The logic flow is depicted in 
Figure i which is a block diagram summary of the procedure. 

Comment 1: As the optimum solution is approached, the gradient 

^ must tend to zero. If this slope is below some acceptable level, then 

we have obtained the solution. If this slope is not below an acceptable level, 
then we have to repeat the procedure. 

Problem Formulation 

Consider a system defined by the n state variables: 


x (t) = 

(n x 1) 


x t (t) 

x 2 (t) 


x (t) 
n 


U) 


These state variables are subject to the m control variables: 


« (t) - 

(m x 1) 


a t (t) 
ot 2 (t) 



(t) 


(2) 


The optimization problem consists of determining the control matrix 
a (t) in the interval t Q to T so as to maximize the function: 
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MAKE INITIAL ESTIMATE 
OF a It) : a' (t) 



USf a' (t) AND INITIAL STATE 
VARIABLE CONDITIONS TO 
CALCULATE NUMERICALLY x' <t) 



PERTURB a' (t) BY AN 
AMOUNT 5 It) TO MEET 
REQUIREMENT OF 
MAXIMUM <t> 


BUTar (t) AND x' (t) 
DO NOT NECESSARILY 
PROVIDE THE DESIRED 
SOLUTION 



5(t) PERTURBATION 
CAUSES PERTURBATION 
IN x (t), 8x (t) 



d <t> : CHANGE IN 0 FOR 
SMALL PERTURBATIONS 
5«<t) ABOUTO* It) 


FOR STEEPEST ASCENT: 
5a(t) THAT MAXIMIZES cty 
FOR GIVEN VALUES OF 
IdP) 2 IS THE DESIRED 


dP : CHOSEN SUCH THAT 
PERTURBATIONS 5 alt) ARE SMALL 
ENOUGH FOR LINEARIZATION 
LEADING TO BOUNDARY 
CONDITIONS FOR ADJOINT 
EQUATIONS TO BE VALID 
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<(> = <t> [Xj (T), x 2 (T), .... 


(3) 


x n (T), T] 


while satisfying the p terminal constraint functions: 


* 

(P x 1) 


♦■llMT), X 2 (T) x n (T) , T] 

^ 2 [Xi(T), x 2 (T) , . . . , x n (T),T] 


* [xj(T) , x 2 (T) 

tr 


x n (T),T] 


0 


(4) 


and the single stopping condition: 


S2 - [x 1 (T),x 2 (T>...,X n (T),T] = 0 


(5) 


and the n system equations: 


f 

(nx 1) 


dxi 

j 

dt 


dx 2 


dt 


dx 


n 


dt 



m 


m 


. ( 6 ) 


fy [ Xj ( t) , X 2 ( t) , . . . , X n (t) , Q! j ( t) , OL 2 ( t) , • > ■ , ( t) , t] 


% 


The optimization problem thus formulated can be solved by using the 
steepest ascent method developed by Bryson and Denham [3] . It is assumed 
that the initial conditions x(t 0 ) are specified. 
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To start the procedure an initial estimate must be made of the control 
variables in the interval between t 0 and the time at which equation (5) is 

satisfied. Denote this estimate by a ^(t) and use it and the initial conditions 

and equations (6) to numerically determine (T) ; i.e., substitute of ^ (t) 

in equations (6) and numerically integrate the resulting equations from the 
initial state until the stopping condition of equation ( 5) is reached . This 

yields the state variables x^ ^ (T) and the final time T resulting from the 

estimate (t) . The values of x^ (T) and T will not necessarily satisfy i' 

[equations (4)] or maximize <p [equation (3)]. Thus, a ^ (t) must be changed 

by an amount 6a to meet these requirements. 

Linear Differential Perturbation Equations 

Towards this objective, consider small perturbations 6a about the initial 
estimates of the control variables. This would lead to a second approximation 
given by: 

a^ (t) = (t) + 6a. , (7) 

where i = 1, 2,...,m. The perturbations cause perturbations in the state 
variables: 

x (2) (t) = x (1) + 6x. (8) 

J J J 

where j = 1, 2, . . .,n. Next, take the variations of equations (6) by first sub- 
stituting equations (7) and (8) into equation (6) and then expanding the right 

side of equation (6) in Taylor's series atiout the nominal estimates a.^ ^ (t) 

and xf ^ (t) : 

J 
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f. |xj 2 ^(t) ,X 2 2 ^ (t) x^ (t) , a i 2 ^ (t) , (t) , . . a^ 2 * (t),t 


.,[■ 

■-[* 


f. |xj 1} (t) , x 2 ( i} (t) x* (t) , a^ (t) , a[ i] (t) , . . . , a^ (t) , t 


n 

z 

H 


8f. 


^ (*> -f '<«] 


m 

- z 

k=i 


af 

5T (“f"* -f ,(t) 

k 


( 9 ) 


where i = 1, 2, . . . , n. The system of equations ( 9) can be written as: 


f ," J’ 


6f ; = f. | xj ^ (t) + 6x t , x 2 ^ + 6x 2 , . . . , x^ ^ (t) + 6x m> a l ^ (!) + i» 


a: 2 ^ (t) + 6a 2 , . . . , a^ ^ (t) + 6a 

1 ' *■ m m 


x^tt), a ( 1 1) (t), a 2 (1) (t) 


i - ■■[ 


xS^(t), x 2 (1) ( t) , 




,2) (t) -x (1) (t, 


3}* Lj=i : 


(2) (t) -^ (1) (t) 


n 

"af. 


m 

af. 

z 

j=l 

1 

ax. 

3 

6x^ (t) 

+ t 

k=i 

— ^ 6a (t) 
3a. k 

k 


( 10 ) 


where i= 1, 2, . n. Thus, the variations in a ^ (t) and xf^ (t) induce a 
variation in f denoted as 6f.. But from the variational calculus, we can write 

l 

by combining equations (6) and (10) [8] : 

(t) 


/<&. \ ,, r 

6f. = 61— r~)= TT 6x. 
i \ dt / dt 1_ i 


or 


r 

n 

r 9f . 


m 

" af. 

|6x. (t) 

= z 

j=l 

1 CO 

J 

6x^ (t) 

- z 

k=l 

— — 6a (t) 
da, k 

k 


( 11 ) 


( 12 ) 
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where i = 1, 2, . . . , n. These equations can be written in matrix notation as: 


^ [fix (t)_ = F (t) 6x (t) + G (t) to (t) 


£[ 6x < t il = 


(n x 1) 


^ [ ax * (t >] 

[**> <*»] 

• I 

£ |f x n (t >] 


6x (t) = 
(nx i) 


6xj (t) 

<5x 2 (t) 


fix (t) 
. n 


to j (t) 
to 2 (t) 


6a (t) = ' 

(m x 1) 


6a (t) 
m 
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F(t) = 
(n x n) 


d Jl 

9£l 

9 f i 


9x 2 

• • • ■ 

9x 


n 

Qh 

9f 2 

9f 2 

9Xj 

• 

9x 2 

9x 

n 

• 

• 

8f 

9f 

8f 

n 

n 

n 

9Xj 

9x 2 

••• 9x 

n 







9ii 

8 f 1 

9 or t 

da 2 

* ” 9 or 


m 

9fi> 


9f 2 

9a t 

9 Qf 2 

... 9a 


G(t) = 
(n x m) 


af af af 

n n . . . n 

da t da 2 ace 

1 i m 


( 17 ) 


( 18 ) 


The [ ]* indicates that the partial derivatives are evaluated for the 
nominal estimates. The system of linear differential perturbation equations 
for 6 x 5 as given by equations ( 13) , play an important role in determining the 
changes in $ , "5, and £2 as caused by control variable perturbations 6a (t) . 
The following section ascertains this role. 

Changes in *, and s Resulting From 
Control Variable Perturbations 

We now wish to determine the total changes d 0. d^ , and d£2 in <P, ^ , 
and £2 , respectively, for small perturbations So: in the control variables about 
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the nominal estimates. To do this we introduce the linear differential equations 
adjoint to equations (13) and defined as [9] : 


[y‘U - -r T royt) 

(i9) 

[V«] = -F T (t) ^ (t) 

(20) 

[V‘] - (t) 

(21) 


where 


1 M ■ 

(nx 1) 


dA 


& 


dt 

dA 


£2 


dt 


dA 


0 n 


L dt J 


( 22 ) 


A (t) = 

0 

(nx 1) 


01 

0 2 


0 n 


(23) 
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(27) 


X fi (t) = 
(n x 1) 


fii 


fi2 


fin 


By definition, the boundary conditions for these equations are given by 
the following: 


V T) 




where 


90 

9x 

(lx n) 


9j£ 

9xj 


90 

_ ■ • • 

9X2 


9 0 

9x 


n 


JL£_ 

9x 

(px n) 


9* t 90j 

9x t 9 x 2 

9^2 9J_i 
9x t 9 x 2 


EL 

9xj 




n 




% 

(28) 


(29) 


(30) 


(31) 


( 32 ) 
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an 

8x 


(33) 


an a_fi_ 

Dxj 9 x 2 " ’ ’ ax n 


The () * indicates that the partial derivatives are evaluated for the 
nominal estimates. If we now take the transpose of equation ( 19) , post- 

multiply by 6x and pre-multiply equation (13) by A T and add the results, the 

following is obtained: 

5?F X <«] + i [}* ^ <*> ° <‘> ^(t) . (34) 

( ixn) (nxl) (lxn) (nxl) ( lxn) (nxm) (mxl) 


We can also write : 


d 



(t) <5x (t) 
dt 


(t) w [ & (t) ] + [% T (t) ] 6x (t> • 


(35) 


If we now substitute equation (35) into equation (34) , the following 
equation is obtained: 



(t) 6x (t) 


;) 


dt 


C (t) G (t) fice(t) 
<P 


(36) 


By integrating equations (36) from t 0 to T, we obtain: 

T r -j T _ 

f dxj(t)fix(t) = / XT (t) G (t) 6a (t) dt (37) 

to L 9 -> t 0 * 


i 


\ 

k 


13 




i 


or 


- T — 


- T — 

A (t) 5x (t) 
0 

t= T 

A^ (t) <5x (t) 


T _ _ 

= f xj (t) G (t) Sc* (t) dt 

f 0 

*0 


( 38 ) 


In an analogous manner, 


r - t — l 

r- t — i 

A (t) Sx (t) — 

1 A^ (t) 5x (t) 

L x = T 


t — tn 


T _ 

/ A^ (t) G <t) 6a (t) dt 
to 

- [ 


(39) 


/ 


xj'(t) 6x(t) 


t = T 


xju) <5 x (t) 


t = t n 


r - T — — 

= J \ (t) G (t) 6a (t) dt 


(40) 


We now consider the functions and £2 that are given by equations 
(3) , (4) , and (5) . The total differentials of these expressions are: 


d d> = 
(lx i) 




dx 


dx + 


9 (ft 
9T 


dT 


d* = dx + 
(pxl) 9x 9T 


dT 


dfi = 

(lx 1) 


M £ + dT 

9x 9T 


* 


(41) 
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where 


dx = 
(nx i) 


But equations (42) can be expressed as (see Reference 9 for a dis- 
cussion of these equations) : 


6x t (T) 
Sx 2 (T) 


+ fdT 


6x (T) 
n 


dx = 6x (T) + f dT 


By substituting equations (44) into equations (41), 


d(p = 6$ + <£ dT 
d^ = 6* + ^ dT 


dfi = 6f2 + ft dT 


where 


6<p = 6x (T) 
(1x1) dx 

6x (T) 


(px 1) 
6fi 

(lxl) 


3 fl — yrriv 

6n = aT" 6x (T) 


( 46 ) 


* = 1 ± + |£ f 


(ixl) 


8t 3x /t = T 


(px 1) \ 0t 9X j t=T 

• = /8fi + an f -\ 

(lx 1) \ at 9X /t = T 


(47) 


If we compare equations (46) with equations (28) , (29) , and (30) , 
we see that : 

60 = aJ (t) (T) 

4 > 


6* = ^ (t) 6x (T) 
6fi = (t) & (T) 


(48) 


Substituting equations (48) into equations (45) and rearranging: 


A^(T) 6x (T) = d0 - 0 dT 
A^ [ '(T)6x(T) = d? - * dT 
aJ (T) 6x (T) = dfi - fi dT 


(49) 
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We can now substitute equations (49) into equations (38) , (39) , and 
(40) to obtain: 



(50) 


These are the desired changes in (p, V , and £2 for the control variable 
perturbations 6a (t) . 


Determination of the Control Variable Perturbations 


We now wish to find the 6a that maximizes d<£ in the first expression of 
equation (50) for given values of: 1 

T 

1) (dp) 2 = J 6a (t) w 6a (t) dt (51) 

to 

2 > « 

3) dfi 


t First, we solve the third equation of equation (50) for dT and eliminate 

it from the first and second expressions of equation (50): 


1. The matrix wis an m x m symmetric weighting matrix chosen to improve 

convergence of the procedure . 


I 

t 0 


xj(t) G (t)6a (t) dt + xj (to) 6x (t 0 ) 


(52) 
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Substituting equation (52) into the first and second expressions of 
equation (50), 


or 


d<^> = J xj (t) G (t) do; (t) dt + xj (t 0 )6x (t 0 ) 


4 > 


+<P 


~ -r ^ ^ (t) G (t) 6a (t) dt + X^ (t 0 ) 6x (t 0 )^ 


d* = J X^(t) G (t) 6a (t) dt + X^ A (t 0 ) 6x (t 0 ) 
to 

+ i ^ (t) G (t) 6a (t) dt + xj (t 0 )6x (t 0 )^ 


:T„ - 


( 53 ) 


(54) 


dcp = f X T (t) G (t) 6a (t) dt + (t 0 ) 6x (t 0 ) 

/ 4>1 <P 1 

l o 


r T - 


(55) 


d* = J (t) G (t) 6a (t) dt + xj t (t 0 ) 6x (t 0 ) 


*1 


(56) 


where 


- T 

X^ (t) 

- T 

= X (t) 

_ ± 

- T 

X (t) 

4> 1 

0 


S2 

(lxn) 

(lxn) 


(lxn) 

- T , 

— T , 


" T , 

x *t (t) 


“ I 

(t » 

(p xn) 

(p xn) 




(57) 
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As stated before, the problem is to find the control variable perturbation 

to maximize for given values of d^ and dp, with di2 = 0. This is accom- 
plished by first forming a linear combination of equation (55) with equations 
( 51) and ( 56) through the use of the Lagrange multipliers v and n. The second 
variation of this linear combination is then formed. It is shown in Reference 
3 that this consideration leads_to the following expression for the desired con- 
trol variable perturbation for 6a (t) : 

*Vl (t)v ] ’ (58) 

(nxp) ( pxl) 


ta <t) = 2 7 w_1 G T (t)L^ 1 (t) 

(mxl) (mxm) (mxn) (n xl) 


where m is a constant and v is a p x 1 row matrix of constants. The expression 
for to (t) can now be substituted into the equations for d* and (dp) 2 [equations 
( 51) and ( 56)] to obtain two equations that can be solved for n and v. If we 
first substitute equation ( 58) into equation ( 56) , we obtain the following after 
some manipulation: 


v = I^ 1 (Ii " 2 M dtf) (59) 

where 

T 

= J 

(pxl) t„ 


dt (61) 

d¥ = d¥ - (t 0 ) <5x (t 0 ) . (62) 

(pxl) 

From equations (58) and (59) , we can write: 


i 2 = / 

( DXP) tn 


A* (t) G (t)w-'G T 


<t) A #1 (t) 


IX « (t) 


G (t) w 


-g t 


(tj A. <t) 

01 


dt 


(60) 
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(63) 


— T 1 

6(X <*> = 2m 

v T = Ii T (I 2 _1 ) T " 2 m (H' T ( I 2 -1 ) T (64) 

where 

d4' T = d* T - 6x T (t 0 ) A^ (t 0 ) . (65) 


xj (t) G (t) w 
0 1 


-1 _ 


r T (t) G <t) w 


If we now substitute equations (58) and (63) into equation (51), the 
following result is obtained: 


(dp) 2 





^ T _ 

+ d^ I 2 1 d^ 


( 66 ) 


where 


i» = / 

(lxl) t 0 


- T — — -1 — T 

A (t) G (t) w G 

01 


(t) A 


01 


dt 


Solving equation (66) for m: 


2m = ± 


N 


h - !i T k - 1 II 

o ^ X — i ^ 

(dp ) 2 - d* I 2 d* 


(67) 


( 68 ) 


Equations (59) and (68) give the solutions for m and v. These can be 
substituted into equation ( 58) to obtain the following: 


- -l- T 


6a (t) = ± w ‘G * (t)f A^ 1 I 2 " 1 h 


+ w -1 G T (t) A T I 2 1 d4> 
'J'l 




(dp) 2 - d4^ T I 2 1 d4> 

Ij - Il T I 2 _1 Ii 


(69) 
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This is the desired control variable perturbation that maximizes d 6 
in the first expression of equation (50) for a given value of (dp)* [equation (51)] , 
given p values of d't [the second expression of equation (50)] and d£2 = 0 in the 
third expression of equation (50). Thus, we substitute equation (69) in the first 
expression of equation (50) or equation (55) to obtain the predicted change in 
<f> for the change in the control variables: 



+ ij I 2 -1 d* + A ^ (t 0 ) 6x (t 0 ) . (70) 


For d^ = 0 and 6x (t 0 ) = 0, equation (70) becomes: 


(71) 


or 

(72) 

As the optimum solution is approached and the terminal constraints are 
met, this gradient must tend to zero. The + sign is used if $ is to be maxi- 
mized and the - sign is used if 0 is to be minimized. 

The control variable perturbations as given by equations ( 69) are now 
added to the initial or previous control variable estimates to yield the new 
estimates: 
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a (t) = a (t) + 6a ( t) 
n p 


( 73 ) 


The new estimates a; (t) are now used in equation (6} 

is repeated until the terminal constraints are satisfied and ^ 

dp 


and the process 
— 0. 


APPLICATION OF THE STEEPEST ASCENT METHOD TO AN 
APOLLO THREE-DIMENSIONAL REENTRY PROBLEM 


Problem Formulation 


The specific problem investigated herein has been studied extensively 
by Colunga using the indirect method termed the "modified sweep method" 

[7] . Basically, the optimization problem consists of determining the roll 
angle program /3 (t) in the interval t 0 to t f which can be used to control an 

Apollo spacecraft to minimize the following function: 


to 

This form is of the classical Langrange or Bolza problem and can be trans- 
formed to the Mayer problem [9] . This is done by introducing one additional 
state variable and one additional differential equation. If this state variable 
is denoted as q (t) , then: 


(L 2 _+ D 2 ) * ~ a 


rsj 

m 


+ A 0 p 


dt 


(74) 


q (t) = [integrand of the function I] 


(75) 


or 


dq 

dt 


( Lr + D ) 


, 2\2 


m 


Ao 


p 2 v 3 


(76) 


22 



where the first term measures acceleration caused by aerodynamic forces 
and the second term measures convective heating experienced by the space- 
craft. We also have: 


q (t 0 ) = 0 


(77) 


From equation (75) we can rewrite equation (74) as: 




4 


q(t) 



(L 2 + D 2 )^ 
m 


+ \ 0 p 2 v 3 


dt 


or 


(78) 


I = q (t) 


(79) 


which is of the form of the Mayer problem. Thus, the system is defined by 
the following state variables ( n = 7) : 


x (t) = 
( 7 x 1) 


h (t) 

e (t) 

A (t) 
v (t) 

y (t) 
A (t) 
q (t) 


(80) 


These state variables are subject to the control variable (m = 1) : 


a (t) =[/? (t)] . (81) 

(lx i) 
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The problem is now to determine a (t) in the interval t 0 ^ t < t^ to minimize 


4> = l 


(82) 


This minimization is to be accomplished while satisfying either of the following 
two cases involving the terminal constraints and the stopping conditions: 

Case 1 — Terminal Constraints (P = 5) : 


_ l 

* = 


l 

*1 

l 

*2 

1 

*3 

1 

*4 

1 

*5 


h(t f ) -h f 
0(t f ) - 9 f 
A(t f )- A f 

y (t f )- y f 
A(t f ) - A f 


(83) 


Single stopping condition is as follows: 


ft = v (t^) - 261 m/ sec ( 856 ft/ sec) 


(84) 


Case 2 — Terminal Constraints (P = 5) : 



11 





*1 


e (t f ) 

- £ f 


11 

*2 


A (t f ) 

” A f 

_ 11 

11 


v (t f ) 

■ ; f 

* = 

*3 

— 

y (t f ) 


11 




*4 


A (t f ) 

- A f 


11 



_ 


*5 




( 85 ) 
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Single stopping condition is as follows: 

_ ll 

fl = h (t f ) - 23 014 m (75 504 ft) = 0 (86) 

where 

hj = 23 014 m (75 504 ft) 

= 24. 1 deg 
A^ = -0.6 deg 
= -44. 3 deg 

= 261 m/sec ( 856 ft/ sec) 

= -29. 4 deg 



The time at which the stopping condition is satisfied is denoted as t^. The 

minimization is also to be accomplished subject to the following differential 
equations of motion constraints: 


dh/dt 


v sin y 

d0/dt 


v cos y cos A/[ (R + h) cos A ] 

dA/dt 


v cos y sin A/ (R + h) 

dv/dt 


G sin y - D 

dy/dt 


[(G cos y)/v] + [v cos y /(R + h)l + (Lcos /3 /v) 

dA/dt 


[-vcosy cos A tan A/(R+h)] - {[L sin/3/(v cosy)]} 

dq/dt 


(L 2 + D 2 ) 2 ~ i , 

+ X 0 P 2 v 3 



m 


( 88 ) 


where 



( 89 ) 


G 

D 

L 

L 

D 

P 


-p*/(R + h) 2 ^ 
p S v 2 C d /(2 m) 

p S v 2 C /(2 m) 

f^C L 

^ p S V 2 C D 


The following initial conditions for the state variables are assumed at t 0 = 0: 


h (t 0 ) 


121 920 m (400 000 ft) 

e (t 0 ) 


0 deg 

A (t 0 ) 


0 deg 

V (to) 

= 

10 668 m/sec (35 000 ft/sec) 

x(t 0 ) 


-6. 5 deg 

A (t 0 ) 


0 deg 

q (t 0 ) 


0 deg _ 


Perturbation and Adjoint Equations 

The equations analogous to equations ( 13) governing the behavior of 
perturbations in the system are given by: 

(j — 

— [<5x(t)] = F(t)6x(t) + G (t)<5of (t) (91) 


where 
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4r [5 x(t)] = 
( 7 x 1) 


dT [6h(t)] 


d 

dt 

d 

dt 

d 

dt 

d 

dt 

d 

dt 

d 

dt 


[<50 (t)] 

[<5A(t)] 

[6v(t)] 

[6y(t)] 

[6A(t)] 

[fiq(t)] 


(92) 


6x (t) = 
(7 x 1) 


6h (t) 
<50 (t) 
<5A (t) 
6v (t) 
6y (t) 
6A (t) 
<5q (t) 


(93) 


5a (t)= [6/5( t)] 
(lxi) 


F (t) = 
(7x7) 


(94) 


af< 

a f i 

afj 

ah 

90 

• • • 

aq 

afj 

af 2 

9f 2_ 

ah 

• 

90 

# 

• 

• • • r\ 

aq 

• 

afj 

af 7 

... af^ 

ah 

90 

aq 


(95) 
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( 96 ) 


The linear differential equations adjoint to equation ( 91) are: 



( 97 ) 


X *ii A *12 

X A • • • 

^21 *22 

• • 

• • 

• • 

A t X 
* 71 *72 



r 
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A_ 


A_ 



£2 1 

• 



a_ 


A_ 

n 2 

- T 

Q 2 

• ' 

= -F (t) 

• 

.. • 


♦ 

• 


• 

• 



A 


A_ 



i2 1 


( 99 ) 


The boundary conditions for equations ( 97) , ( 98) , and ( 99) are given by: 


( T > 

<$> 

<7x 1) 


9£ 

ah 

9 0 

a 9 



T 


(100) 


V T > 

( 7 x 5) 



9*j 

9* 2 

a*s 

* • • * 

ah 

9h 

ah 

9*j 

9* 2 

. . . a* s 

80 

90 

9 0 

♦ 

• 

• 

• 

• 

• 

« 

• 

• 

a*j 

9*2 

9* 5 

8q 

9q 

• • • p. 

9q 


( 101 ) 
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For Case 2 (CII) : 


[V T> ]cn = 

(7x 1) 


0 0 0 0 0 
1 0 0 0 0 
0 10 0 0 
0 0 10 0 
0 0 0 1 0 
0 0 0 0 1 
0 0 0 0 0 _ 


1 
0 
0 
0 
0 
0 
0 

Total Differentials of and n 



(7x1) 



(7 X 5) 


0 

0 

0 

0 

0 

0 

1 


(106) 


(107) 


(108) 


The desired changes inj>, and 12 for the control variable per- 
turbations are given by (with 6x (t 0 ) = 0) : 


d <p 



(t) G(t) 6a (t) dt + <j) 


dt 


(109) 
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(110) 


T 

dit = f (t) G (t) 6a (t) dt + ^ dt 

4 

dS2 = f X J (t) G (t) 6a (t) dt + si dt . (Ill) 


If we solve equation (111) for dt and substitute the result in equations (109) 
and ( 110) , 


_ rn 

d<f> = f X , (t) G (t) 6a (t) dt 

4 


d* = / X* (t) G (t) 6a (t) dt 

t * X 


where 


- T 

x" (t) 

<pi 


- T 

<P 


$ X I (t) 

a n 


- T 

V <*» 


- T 

X (t) 

* 


* - T 

l (t> 


* M + a* -A 

y \at 8x It 

(lXl) _ V 

; Id* d* 

" \8t + 8x f /t = 

(5x1) _ 

fi + 8« f ) 

(lx if* 8X /t = 


= T 


(112) 

(113) 


(114) 


(115) 
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For Cases I and II, we have : 




(116) 



( 117 ) 


( 118 ) 




( 119 ) 
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Control Variable Perturbations 


The 6a_(t) that minimizes dcp in equation ( 112) subject to the terminal 
constraints d4s the single stopping condition dft = 0, and the constraint (dp) 2 
on the magnitude of the control change is given by: 


6a (t) = - w - 1 G T (t) (x 0l - I 2 _1 X t ) N 


(dp) 2 - d>fr T l 2 x d^ 

Is ~ Ii T h ' 1 Ii 


+ w" 1 G T (t) X t I 2 -1 d* 


(120) 


where 


T ' 

/ * 

o L 


II = / (t) G (t) w -1 G T (t) X (t) 

(5x1) n 9 


dt 


Is = 

(5x5) 0 


T r 

I VW a(t)w-‘o T (t) x (t) 

o L 


dt 


( 121 ) 


h= J I X* (t) G(t) w -1 G x (t) \ (t) 

(1x1) o 


dt 


We can now substitute equation ( 120) into equation ( 112) to obtain the change 
in (p for the change in the control variable: 


d<p = - ^ [(dp) 2 - d* T I 2 _1 d* j [l 3 - I 2 _I Ijl + Ii T I 2 _1 d* . ( 122) 
For dF = 0, this becomes: 
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(123) 


g-n-«. T 


ir 1 ii 


We now add 6 qj (t) to the initial or previous estimate to yield: 


o; N (t) = “p (t) + da (t) . (124) 

The new estimates a^ (t) are now used in equation ( 88) and the process is 

repeated until terminal constraints are met and 4^- — ■ 0. 

dp 


STUDY RESULTS 
Computer Program Development 

The study of the application problem herein has resulted in the 
development of a highly flexible computer program that can be adapted to 
other trajectory optimization problems. The Marshall Vehicle Engineering 
Simulation System (MARVES) programming system is used extensively in 
this development. The use of the MARVES programming language provides 
the added flexibility of specifying program statements directly related to 
the application problem. It also provides an easy means of modifying the 
developed program. A complete description of the program is given in the 
appendix. 


Application Results 

In an application problem of this nature, several parameters can 
be studied to establish their relative importance and/or effects on the 
minimizing or maximizing function. The Case 1 problem involving the 
five terminal constraints and the velocity stopping condition was selected 
for various parameter perturbation runs. These various runs are summar- 
ized in Table 1. 

The standard deviation of the corrections dec (t) to the control 
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TABLE i. PARAMETER PERTURBATION RUNS FOR THE CASE I APOLLO APPLICATION PROBLEM 


** 1 > 

50.72 

8 . 76 
63.36 

13.99 

4.87 

21.24 

6.23 

10.87 

4.41 

5.23 

4 . 23 
5.98 
12.26 
15.50 
2.64 
3.09 

2.59 

2.94 

2.15 

2.58 

3.38 

5.38 
10.00 

9. 19 

7.38 
3.40 
11.97 
5.80 
6.44 
4.54 

cx 

T 3 

0.01948 

0.01019 

0.02777 

0.00477 

0.0002 

0.0045 

0.0009 

0.0253 

0.00012 

0.00021 

0.00014 

0.00026 

0.00082 

0.0012 

0.00004 

0.00005 

90000 ’0 
50000*0 

0.000035 

0.000049 

0.000077 

0.000164 

0.000478 

0.000393 

0.000322 

0.000098 

0.00072 

0.00021 

0.00034 

0. 00015 

TJ 

i a 

T > 

- 0.0526 

- 0. 1398 
- 0.0823 

ao 

© 

cn 

1 — t 

© 

\ 

- 0. 1398 
- 0 . 1311 

- 0.1398 

- 0.1208 

- 0. 7226 

- 0.7195 

- 0.7468 

- 0.7395 

- 1. 1377 

- 2. 1606 

- 1.362 

- 1.5015 

- 0.1398 

- 0.1382 

- 0.7468 
- 0.7515 
- 0.7684 
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- 1.891 
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- 1 . 760 
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- 2.70 
- 2.22 
- 2.91 


— 

m 

t — 

m 

CD 
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- 0.0172 

1.93 

- 0.0137 

- 0.0034 

- 0.0091 

- 0.0069 

- 0.262 

- 0.0275 

- 0.0356 

- 0.027 
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- 0.045 
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- 0.0061 

- 0.0139 
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- 0.0187 
- 0.0212 
- 0 . 0237 
- 0.0255 
- 0.0266 
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- 0;0282 
- 0.0287 
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- 0.0294 



0.532 

0.3159 

0.5375 

0.253 

0 . 063 
0.083 

0 . 126 
0.189 

- 0.0650 

- 0.069 

- 0.068 

- 0.071 
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- 0.0721 

- 0.065 

- 0.068 

0.0316 

0.0365 
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- 0.0360 
- 0.0361 


co 

00 

© 
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0.0009 

0.0017 
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0.00006 

0,00002 

0.00006 
0.000005 
- 0. 000075 
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60000*0 

0.000030 
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- 0. 000054 
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0.000082 


CM 

© 

© 

d 

i 

- 0.0146 

0.0077 

- 0.0116 
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- 0.0015 
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parameter is computed at the end of each cycle for a given run. This pro- 
vides a quantitative measure of the overall magnitude change in the corrections 
from one cycle to the next. This quantity is denoted in Table 1 as cr^ . 

It is easily seen from Table 1 that the steepest ascent method is 
sensitive to the selection of the parameters k 0 , dp 0 , and kj. These results 
indicate the need for a more rigorous selection criteria for these parameters. 
It should be pointed out that thek 0 , dpo, and kj values used are not necessarily 
the values to use for other application problems. 


APPENDIX A. COMPUTATIONAL SOLUTION PROCEDURE 
FOR THE STEEPEST ASCENT METHOD 


This appendix lists the various computational steps in the implementa- 
tion of the steepest ascent method. These are the specific steps to be taken 
to obtain the minimizing solution for the Apollo application problem discussed 
earlier. A computer program listing containing the computational details is 
also given. Typical output results for the application problem are included. 

It should be pointed out that the program listing is for the deck operation 
on the IBM 7094 computer. A similar deck exists for use on the UNIVAC 1108 
computer. 

It is noted that both programs have double precision capability. The 
SC-4020 plotting procedures are also used extensively in each program. 

Step 1 

Integrate equation (88), given initial estimates [Qf (t)l of control 

p 

variable and initial conditions [x(t 0 = 0)] , until stopping condition of equation 
(84) or equation (86) is satisfied. The time at which this occurs is denoted 
as T. Store the state variable values between t = 0 and t = T. 

Step 2 

Integrate equations ( 97) , ( 98) , and ( 99) backwards from t = T to 

determine the adjoint variables X^ ( 7 x 1) , X^ (7x5), and ( 7 x 1) . 

The matrix F (t) is evaluated on the nominal path by reference to the stored 
values of the state obtained in Step 1. Thus, F (t) is a time varying array 
of coefficients. 


Step 3 


Calculate X 

-T - * 

and X G 
^ 1 


1 ’ 


Xy A as given by equation 


114) . Then form X ^ G 
<t> 1 
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Step 4 

Carry out backwards integrations to obtain I 1# I 2 , and I 3 [equation ( 121)] . 

Step 5 

Print out the values of q, ^5 achieved by the nominal 

trajectory. These values result from using the initial control estimates 

otp (t) . This is the first solution. 

Step 6 

Select dpo so as to obtain a reasonable value of dpjj/T, a mean square 
deviation of the control from the nominal to the next step. 

Step 7 

Select the changes d^ 1( dtf 2 , . . . , d^ 5 so that 



t = T 

within the limitation that: 

dp = ^ dpo - d'J' ^ I 2 1 d^J s 0 

and where ^ refers to Case 1 [equation (83)] or Case 2 [equation (85)] for 
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rsj rss 

the terminal constraints. If the 4' *s make dp negative, reduce d4> by a 
constant factor so that dp = 0 (or nearly vanishes) ; i.e. , the limitation on 
dp becomes: 

dp = kj) ^dpo - kf d^ T I 2 -1 dtyj > 0 

where 0 < k] < 1. It is noted that two cases can occur; namely, dp > 0 or 
dp < 0. If dp > 0, we would proceed to Step 8. However, if dp < 0, then 
k t is selected such that dp > 0 and then control is transferred to Step 8. 

Step 8 

Use equation (120) to determine 6a (t) and add it to a (t) to form 

P 

a (t) + 6a (t) for the next approximation. 


Step 9 




Examine the 4* l5 4^, • • • > 4' 5 values from Step 5 to see how close to 
zero they are. Examine [equation (123)] to see if the solution is sufficiently 
close to the minimizing solution. 


Step 10 

For 4» s* 0 satisfied, if ^ is not sufficiently close to zero, return 

to Step 1 and use a (t) + 6a (t) and repeat the computational cycle. 

P 

Step 11 


d (h 

Terminate when 4> = 0 and ^ — 0 

dp 


Steps 1 through 5 yield all the information for the first cycle; i.e. , we 
have the optimum payoff function q (t) which results from using the control 
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variable estimate o' (t) 
P 


The information required in the following cycles is 
given in Steps 6 through 11. A value for the parameter k 0 is selected to obtain 


values for d'k which are closer to zero. A constant k 0 is maintained in all 
cycles from the second value on. In checking ^ , if oscillations occur, then 
dpo is reduced by an order of magnitude; i.e., dp = 0. 01, 0. 001, 0. 0001, etc. 


Example of Plot Output Results for Case I Application 
Problem, Run Number 6 (Cycle 2) Data 


H: 

Altitude (Fig. A-l) 

Theta: 

Longitude (Fig. A -2) 

Delta: 

Latitude (Fig. A -3) 

V: 

Velocity (Fig. A-4) 

Gamma : 

Angle of Attack (Fig. A -5) 

A: 

Heading Angle (Fig. A -6) 

F7A: 

Acceleration Component of Q (Fig. A-7) 

F7H: 

Heating Component of Q (Fig. A-8) 

Beta: 

Roll Angle (Fig. A-9) 



100 


200 


300 


400 


Figure A -2 


TIME 


. Longitude versus time. 
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DELTA 



TIME 


Figure A -3. Latitude versus time. 
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TIME 


Figure A-5. Angle -of -attack versus time. 
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F7A 



Figure A-7. Acceleration component versus time. 
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HU 



Figure A -8. Heating component versus time. 
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Computer Listing of Steepest Ascent Optimization Program 
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2 1 J TeMP(IJKL + 7)aTLPHI(IJKL)-TE.lP(IJKL) 

CALL MAT MUL(BTEMP.T6,TEr4P<8>. -1,7,1 
ATAbLE ( JJ ) *.< * BTEMPt i ) « OPC 

CALL BATMUL(ATEMP.AIV2.PSIrt,S,5,l ) 
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END 
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Typical Output Results for Case I Application Problem, 
Run Number 8 (Cycle 1) 
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